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ABSTRACT 


Three  different  algorithms — the  power  series,  the  asymp¬ 
totic  series,  and  the  recurrence  relation  method — are  investi¬ 
gated  with  special  attention  to  the  single  (10“^“^)  and  double 
precision  (10  of  Computer  Data  Corporation  (CDC)  6000 
Series  computers.  The  final  accuracy  of  each  method  depends 
partly  on  the  magnitudes  of  the  largest  and  smallest  terms  when 
floating  point  additions  are  involved.  Another  consideration 
is  the  number  of  terms  required  for  each  algorithm.  Combina¬ 
tion  of  all  considerations  leads  to  a  partitioning  of  the  order- 
argument  domain  into  partially  overlapping  areas  in  which  each 
algorithm  is  most  appropriate.  A  wedged  area  not  covered  by 
any  of  the  algorithms  remains  for  large  order  and  argument 
of  approximately  equal  size. 

Orders  and  arguments  up  to  1024  were  investigated  and 
checked  where  possible.  A  FORTRAN  IV  program  in  the  form  of 
an  external  function  is  included. 


ADMINISTRATIVE  INFORMATION 

The  particular  work  addressed  herein  was  supported  by  the  Naval 
Ship  Systems  Command  (037)  under  Program  Element  25684,  Project  S4628, 
Task  Area  S4628-019,  Naval  Ship  Research  and  Development  Center  Work 
Unit  1932-010. 


INTRODUCTION 

Three  well-known  algorithms  to  determine  the  Bessel  function  (X) 

J-i 


are: 


(1)  The  power  series  expansion^  for  X/L  not  too  large. 


Abramowitz,  M.  and  I. A.  Stegun  (Ed.),  "Handbook  of  Mathematical  Functions," 
National  Bureau  of  Standards,  U.S.  Government  Printing  Office,  Washington, 
D.C.  (1954),  p.  360. 


2.  The  as5nnptotic  series  expansion  for  X/L  large. 

3  4 

3.  The  recurrence  relation.  ^ 

The  power  series  converges  always.  However,  for  large  values  of 
X/L,  a  large  number  of  terms  is  necessary  to  reach  a  given  precision.  More 
seriously,  the  terms  initially  increase  in  magnitude,  starting  to  decrease 
only  after  going  through  a  maximum.  When  the  terms  are  summed  by  employing 
floating  point  techniques,  using  a  fixed  number  of  significant  digits, 
the  accuracy  of  the  final  answer  decreases  by  the  same  order  as  the 
order  of  magnitude  of  the  maximum  term.  In  other  words,  the  final 
accuracy  corresponds  to  the  number  of  digits  available,  counted  from  the 
most  significant  digit  of  the  maximum  term. 

The  asymptotic  series  has  an  additional  problem.  Since  it  eventu¬ 
ally  starts  to  diverge,  the  minimum  term  may  not  be  small  enough  to 
reach  the  required  precision.  In  that  case  this  method  cannot  be 
applied  without  further  refinements.  Even  if  the  required  precision 
can  be  reached,  the  asymptotic  series  may  still  be  inappropriate  because 
of  the  number  of  terms  involved  or  because  of  the  magnitude  of  its 
maximum  term,  which  affects  the  accuracy  the  same  way  as  described  for 
the  power  series  in  the  previous  paragraph. 
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Abramowitz,  M.  and  I. A.  Stegun  (Ed.),  "Handbook  of  Mathematical  Functions," 
National  Bureau  of  Standards,  U.S.  Government  Printing  Office,  Washington, 
D.C,  (1954),  p.  364. 
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Abramowitz,  M.  and  I. A.  Stegun,  "Generation  of  Bessel  Functions  on  High 
Speed  Computers, "Mathematical  Tables  and  Aids  to  Computation,"  Vol.  II 
(1957),  pp.  255-257. 
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British  Association  for  the  Advancement  of  Science,  "Mathematical  Tables," 
Vol.  X,  Bessel  Functions,  Part  II,  Functions  of  Positive  Integer  Order, 
Cambridge  University  Press,  Cambridge  (1952),  p.  XIV. 
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The  algorithm  based  on  the  recurrence  relation  appears  to  be  the 
least  sensitive  method.  The  only  disadvantage  is  that  the  normalizing 
factor  of  the  algorithm  consists  of  the  sum  of  all  even  order  Bessel 
functions,  down  from  approximately  twice  the  value  of  the  argument.  For 
large  arguments,  this  number  may  become  prohibitive.  The  effect  of  such 
long  sums  on  the  accuracy  of  the  method  has  not  been  determined  explicitly 
in  this  paper,  but  some  checks^  showed  that  the  absolute  accuracy  was 
always  within  three  decimals  of  the  machine  precision.  Usually  it  was 
much  better. 

The  program  given  in  Appendix  A  applies  bounds  that  are  based  on 
a  number  of  considerations  described  in  detail  in  the  chapters  that 
follow.  These  considerations  concern  the  number  of  terms  required  and 
the  magnitudes  of  the  maximum  and  minimum  terms.  The  precisions  for 
which  the  bounds  are  derived  correspond  to  the  single  (10~^^)  and 
double  (10  precision  modes  of  the  Control  Data  Corporation  (CDC) 

6000  series  computers. 

Repeated  use  will  be  made  of  a  form  of  the  Stirling^  asymptotic 
approximation  to  the  factorial: 


k; 


y2rt 


(1) 


This  approximation  is  already  accurate  to  2  percent  when  K  =  4, 


Hayashi,  K,  "Tafeln  der  Besselschen,  Theta,  Kugel  and  anderer  Funktionen," 
Springer  Verlag,  Berlin,  Germany  (1940). 

Feller,  W. ,  "An  Introduction  to  Probability  Theory  and  Its  Applications," 
John  Wiley  and  Sons,  Inc.,  New  York,  Vol.  I  (1968),  p.  52. 
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POWER  SERIES  EXPANSION 


The  power  series  expansion  for  J  (X)  is 

Li 


Jl(x) 


f  (E  f  E  ^ 

V  2  y  V.  2  J  V  2  y 

l:  ■  i:(i+l):  ■^2:(2+l): 


(-1) 


^x  \ 
K  'v.  2  y 


2K+L 


k:  (k+l) : 


(2a) 


(Jj 

l: 


1- 


^  X  -f 
V  2  y 
l+(l+L,^ 


+  (-l) 


K 


.  ..2(K-1)+L  ,  2 

(I;  kil 

(K-l): (K-l+L);  K(K+L) 


(2fa) 


The  desired  precision  q  is  reached  when  the  last  term  T  becomes 

K 

less  than  that  precision.  The  condition  is  thus 


T  ^  q 
K  ^ 


(3) 


where 


T 


K 


(if! 

K.'  (K+L).' 


(!)■ 


2  K+L 


K+K+L 


2rt  ^  (K+L)*^+^ 


(4) 


Stirling's  approximation  to  the  factorial  has  been  used  to  obtain  the 
final  expression  in  Equation  (4). 

Instead  of  solving  for  the  K  that  satisfies  Equation  (3),  it  takes 
in  general  less  computer  time  to  test  if  the  last  term  T  has  indeed 
reached  the  desired  precision  q. 
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The  precision  of  the  last  term  does  not  guarantee  the  accuracy  of 
the  final  answer.  The  magnitude  of  the  maximum  term  may  be  several  orders 
of  magnitude  larger  than  one,  so  that  an  equal  amount  of  accuracy  is  lost 
in  the  final  answer  of  this  alternating  sign  series  if  floating  point  tech¬ 
niques  are  employed.  Therefore,  the  term  with  the  largest  magnitude  must 
be  determined.  This  term  occurs  when  the  incremental  multiplier  in  the 
last  term  in  Equation  (2b)  first  becomes  less  than  one: 


K(K+L) 


This  implies  that 


where 


X^/4  <  +  K  L 


K  >  K, 


Ki  =  -  +  i-  /L®  + 


The  magnitude  T  of  this  maximum  term  follows  from  Equations  (4)  and  (7). 


2K,+L 


TKi  = 


(?r 


C?) 


2Ki+L 


2jt[Ki(Ki+L)]^l'^  ^(Ki+L)^  2x[(-  ^L+  |/l^+X^)  (|l+  |/l®+X^)]^’-'^  ®  (Ki+L)^ 


,2K,+L 


(  X  _ 

K2  J  exp(-lrb/L‘^+X^+L)  ^  exp[Xv/(L/X)^  +  1] 

.  -n2Ki+1  ^ ^ _ 

2rt  J  (^L  +  rtX[L/X  +  /(L/X)^  +  1] 


exp(-lrb/L‘^+X^+L) 
,2K,+1  _ 


The  behavior  of  the  magnitude  of  the  maximum  term  T  of  Equation  (8) 

K-i 


is  shown  in  Table  1. 
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TABLE  1  -  MAGNITUDE  T^  OF  MAXIMUM  TERM  IN  POWER  SERIES  OF  (X) 

JL 


L 

1 

10 

100 

1000 

1 

0.5  (K=l) 

600  (K=5) 

10^^  (K=50) 

10^®^  (K=500) 

10 

10“®  (K=l) 

6  (K=2) 

10^°  (K=45) 

10^3  0  a=495) 

100 

lO"’-®'^  (K=l) 

CD 

00 

0 

10®°  (K=21) 

10^®°  (K=452) 

1000 

0  (K=0) 

0  (K=l) 

10“®®®  (K=2) 

10®®®  (K=207) 

A  more  detailed  table  is  given  in  Appendix  B.  If  a  loss  in  accuracy 
of  no  more  than  two  or  three  decimals  is  required,  the  maximum  term  must 
be  less  than  10^*®  or  approximately  300, 


ASYMPTOTIC  SERIES  EXPANSION 


.  2 


The  asymptotic  series  expansion  for  J  (X)  is 

Li 


JlW  [P(L,X)  cos  Y  -  Q(L,X)  sin  Y] 


(9) 


with 


X  -  (^L  +  i)  It 


P(L,X)  ~  1-  ^  JTCSXp  4J  (8X)^ 


Qa,x) 


•%/ 


(4L^-1^) 

IJ  8X 


(4L^-1^) (4L^-3^) (4L^-5^) 
3;  (8X)3 


Consider  the  term 


(11) 
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R^a^x)  = 


[(2L)^-l^][(2L)^-3^] -  [(2K-1)^  -(2L)‘ 

K.'  (8X)^ 


(12a) 


,  ^  I (2K-1)^  -(2L)^| 


(12b) 


These  terms  decrease  as  long  as 


(2K-1)®  -(2L)' 


This  yields 


4L^  -  4K®  +  4K  -  1  -  8XK  <  0 


and  4K®  -  4K  +  1  -  4L®  -  8XK  <  0 


Solving  Equation  (14)  gives  for  the  interval  of  K  for  which  the  terms 
decrease  from  maximum  to  minimum  magnitude 


Ka  <  K  <  Ka 


where 


=  -  i-(2X-l)  +  l-V  (2X  -1)^  +  4L®  - 


K3  =  ^(2X+1)  +  |■V(2X+l)®  +  4L^  -  1 


It  is  obvious  that  Kq  <  L  <  K3,  so  that  K3  -  L  and  L  -  Kg  are  both 
positive,  a  fact  that  will  be  used  in  the  four  following  equations. 
From  Equation  (12a)  it  follows  that  the  term,  if  K  >  L,  can  be 


written  as 
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D  /T  v^  -  [ (1+2L) (3+2L) . . (2K-1+2L) ] f (2L-1) (2L-3) . > 3, 1. 1. 3. . (2K-1-2L) ] 
-  K  (8X){^ 


(2K+2L):  Li  (2L)J  (2K-2L); _ 1 

(2L):  2X  (K+L):  2L  l:  2K-L  (K-L):  KJ  (8X)J^ 


(2K+2L);  (2K-2L)J 
(K+L).'  (K-L): 


1  1 
^  K 


^  P  (K  >  L) 


(18) 


Similarly,  for  K  <  L 


Rj^(L,X) 


[  (2L  +  1)  (2L+3) . . (2L+2K-1) ] [ (2L-1) (2L-3) . . (2L-2K+1) ] 
k;  (8x)^ 


(2Lr|-2K)  ]  L;  (2L);  (L-K):  1  1 

(2L).'  2K  (L+K):  (2L-2K)*  2K  Lj  Kj  (8X)1^ 


(2L+2K).'  (L-K)i  1  1 

(L+K).'  (2L-2K):  2^  Kj  W 


(K  <  L) 


(19) 


Using  Stirling's  formula,  we  find  asymptotically 


and 


(T  (K+L)^'^^ 


1 

v/2^K  ^  K2Xe 


(21) 


By  substituting  Kg  and  Kg  from  Equations  (16)  and  (17)  in  Equations 
(21)  and  (20),  respectively,  we  obtain  the  maximum  and  minimum  magnitude 
of  the  asymptotic  terms.  The  minimum  term  should  be  within  the  precision 
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required  of  the  answer,  while  the  accuracy  is  determined  by  the  number  of 
decimals  lost  in  floating  point  arithmetic,  i.e,,  by  the  order  of  magni¬ 
tude  of  the  maximum  term  R  . 

*^2 

If  either  of  these  requirements  cannot  be  met,  a  different  algorithm 
must  be  used. 

Tables  of  terms  with  maximum  and  minimum  magnitudes  are  given  in 
Appendix  C. 


RECURRENCE  RELATION 

3  4 

The  recurrence  relation  for  Bessel  functions  is  ^ 

-  f 

with  the  normalization  constraint 

00 

Jq(X)  +  2  ^  J2jj(X)  =  1  (23) 

N=1 

The  algorithm  based  on  this  recurrence  relation  starts  with  setting 
Jj^j^(X)  equal  to  zero  and  equal  to  a  (small)  constant.  A  running 

count  is  kept  of  the  normalizing  sum.  Equation  (23).  The  final  answer 
is  obtained  by  dividing  the  value  resulting  from  the  recursive  iteration 
by  the  normalizing  factor.  The  method  is  remarkably  insensitive  to  the 
starting  point  and  the  starting  value.  It  is  advisable  to  start  the 
recursion  at  an  order  L  for  which  J  (X)  equals  approximately  the  desired 

Li 

precision  q,  which  is  usually  the  machine  precision  available. 

This  starting  order  is  readily  estimated  from  the  power  series 
expansion  given  by  Equation  (2b): 
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(24) 


Cf) 


l: 


h=  rt")  = 

9TfT.  ^  2L  y 


72^ 


and,  again  from  Equation  (2b), 


(25) 


The  latter  inequality  ensures  that  the  series  in  parentheses  in  Equation 
(2b)  converges  from  the  start.  Since  the  normalizing  factor  of  Equation 
(23)  has  to  be  carried  all  the  way  down  to  the  zero  order  function  Jq(X), 
a  too  high  starting  order  may  be  computationally  unacceptable;  however, 
in  that  case  an  asymptotic  expansion  may  already  perform  well.  An 
additional  advantage  of  the  recurrence  relation  method  is  that  function 
values  for  all  integer  descending  orders  may  be  obtained  by  the  same 
effort. 

Estimates  based  on  Equation  (24)  for  the  starting  order  L  of  the 
recurrence  relation  of  Equation  (22)  are  given  in  Appendix  D.  From  the 
results  of  Appendix  D  it  follows  that  L=1.4X  +  25  (for  single  precision 
q=10  )  or  L=1.6X  +  40  (for  double  precision  q=10  )  are  good  estimates 

for  starting  orders  of  the  recurrence  relation  algorithm  for  values  of  X 
up  to  200.  In  general,  one  would  not  like  to  extend  the  recursion  over 
more  than  200  terms,  in  the  first  place  because  of  possible  loss  of 
accuracy  involving  the  normalizing  sum  of  Equation  (23),  and  in  the  second 
place  because  of  computational  efficiency.  From  the  requirement  that  L 
be  less  than  200,  we  can  derive  the  requirement  that  X  be  less  than  125 

—  1 4 

(for  single  precision  q=10  ),  or  that  X  be  less  than  100  (for  double 

.  .  ._-39. 

precision  q=10  ). 
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How  accurate  the  obtained  results  are  for  orders  close  to  their 
starting  order  could  not  be  established  exactly,  but  it  is  probably  close 
to  or  better  than  the  machine  precision.  Where  the  values  obtained  by  the 
recurrence  relation  algorithm  could  be  checked  against  Reference  5,  it  was 
generally  found  that  they  were  more  accurate  than  the  corresponding  values 
obtained  by  the  power  series  or  asymptotic  series  algorithms,  which 
always  suffer  from  the  loss  of  accuracy  due  to  the  magnitude  of  the  maxi¬ 
mum  term  occurring  early  in  this  alternating  series,  as  was  explained  in 
the  previous  chapters. 

The  high  accuracy  obtained  by  the  recurrence  method  over  the  full 
range  of  orders  L,  together  with  the  simplicity  of  the  algorithm  itself 
(three  multiplications  and  two  additions  per  step,  following  by  two  memory 
exchanges),  makes  the  recurrence  relation  algorithm  very  attractive  for 
the  computation  of  J  (X).  For  this  reason,  some  of  the  boundaries  indi- 
eating  when  the  power  series  expansion  and  when  the  asymptotic  series 
expansion  can  be  used  may  be  relaxed  with  respect  to  the  results  from 
Appendices  B  and  C.  The  boundaries  actually  used  in  the  program  are  given 
in  the  next  chapter. 

ORDER-ARGUMENT  DOMAIN  COVERED  BY  ALGORITHMS 
The  three  algorithms  into  which  the  calculation  of  the  Bessel  func¬ 
tion  J  (X)  is  divided  occupy  the  following  regions  in  the  order- argument 

Lt 

(L-X)  domain: 

1.  Power  series  expansion  when 

X  less  than  1  or  X  less  than  -l-L  . 
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2.  Asymptotic  series  expansion  when: 

X  greater  than  50  and  X  greater  than  (for  mach.prec.lO 

or  X  greater  than  30  and  X  greater  than  (for  mach.prec.lO 

3.  Recurrence  relation  method  when: 

X  less  than  100  and  L  less  than  200  (for  mach.prec.lO  ) 

-14 

or  X  less  than  125  and  L  less  than  200  (for  mach.prec.lO  ). 

The  coverage  of  the  different  algorithms  is  shown  graphically  in 
Figure  1.  A  wedged  area  uncovered  by  any  of  the  algorithms  remains  for 
large  orders  and  arguments  of  approximately  equal  magnitude. 

Outside  the  rectangular  area  in  the  L-X  domain  for  which  the  recurrence 
relation  algorithm  seems  indicated,  the  boundaries  for  the  power  series  and 
asymptotic  series  expansions  could  be  tightened  again,  but  this  has  not 
been  attempted  in  the  program  presented  in  Appendix  A.  Instead,  the  pro¬ 
gram  notifies  the  user  of  any  entries  into  this  forbidden  zone,  but  returns 
with  a  Stirling  approximation  to  the  first  term  in  the  power  series  -  which 
may  be  completely  wrong  -  as  an  exit  value. 

CHECKING  THE  RESULTS 

To  check  the  accuracy  of  the  algorithms,  Bessel  functiorsJ  (X)  were 
calculated  by  the  recurrence  relation  method  for  arguments  from  1  to  10 
in  steps  of  1,  and  from  10  to  100  in  steps  of  10.  The  region  covered  in 
the  L-X  domain  in  this  way  is  pictured  as  the  hatched  area  in  Figure  2. 

Also  indicated  in  Figure  2  are  the  points  in  the  L-X  domain  for  which 
values  of  the  Bessel  function  were  calculated  by  the  power  series  expansion 
or  by  the  asymptotic  series  expansion.  The  values  of  L  and  X  were  chosen 
close  to  the  boundaries  assumed  by  the  FORTRAN  program  of  Appendix  A. 
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The  required  accuracy  —  three  decimals  less  than  the  machine  pre¬ 
cision — was  indeed  obtained;  this  is  in  accordance  with  the  predictions 
the  previous  chapters.  The  accuracy  obtained  by  the  recurrence  relation 
method  was  generally  better. 


Some  of  the  results  obtained  are  given  below. 


J^(l)  (Ref.  5) 

0.76519 

76865 

57966 

55144 

97175 

26103 

(recurrence  rel.) 

2612 

(power  series) 

2604 

J^(IO)  (Ref.  5) 

-0.24593 

57644 

51348 

33519 

77608 

62485 

(recurrence  rel.) 

6253 

(power  series) 

4585 

Ji(l)  (Ref. 5) 

0.44005 

05857 

44933 

51595 

96822 

03719 

(recurrence  rel.) 

0372 

(power  series) 

0367 

Ji(40)  (Ref.  5) 

0.12603 

83180 

37584 

99920 

56027 

21839 

(recurrence  rel.) 

2185 

(asympt. series) 

2179 

Ji(50)  (Ref.  5) 

-0.09751 

18281 

25175 

13766 

14589 

53873 

(recurrence  rel.) 

53721 

(asympt. series) 

53782 

J^6(50)  (Ref. 5) 

0.00489 

81607 

77813 

78173 

17342 

69265 

(recurrence  rel.) 

68979 

(asympt. series) 

69234 

Jsod)  (Ref.  5) 

•041 

3482 

8697 

(recurrence  rel.) 

3482 

86979 

42514 

(power  series 

3482 

75946 

71184 

•O36 

3241 

50085 

84477 

63106 

(recurrence  rel.) 

3241 

50085 

84477 

63102 

(power  series) 

3241 

82782 

* 

From  a  60-term  power  series  expansion;  also  from  a  recurrence  relation 
calculation  starting  at  J^g^(64). 
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The  results  given  here  show  for  example  that  J^(IO),  when  computed  by 
using  the  power  series  expansion  algorithm,  is  accurate  only  to  25  decimal 
places.  This  is  consistent  with  the  fact  that  the  maximum  terms 
(T4,=T  =678=10^*®^,  see  Table  4  in  Appendix  B)  decrease  the  accuracy  by 
almost  3  decimal  places.  Note  that  the  program  of  Appendix  A  would  have 
chosen  the  recurrence  relation  method  to  compute  J^(IO),  giving  a  result 
that  is  accurate  to  28  decimal  places. 

Taking  as  another  example,  we  can  derive  from  the  chapter  on 

coverage  of  the  L-X  domain  that  the  program  would  again  have  chosen  the 
recurrence  relation  algorithm.  The  asymptotic  series  algorithm  would  in 
this  case  have  been  slightly  more  accurate,  but  both  results  are  good  to 
more  than  26  decimal  places  as  was  desired. 

J„„(l)  begins  to  show  a  difference  at  the  fifth  significant  digit  of 
the  power  series  result.  This  is  because  of  the  fact  that  the  algorithm 
takes  one  more  term  after  detecting  that  the  machine  precision  has  been 

“B 

reached,  which  in  this  case  is  the  very  first  term.  The  third  term  is  3x10 
times  the  first  term  and  this  accounts  for  the  difference. 

The  last  result,  J  (64),  again  would  have  been  calculated  by  using 
the  recurrence  relation  method;  according  to  the  rule  of  thumb  given  in 
Appendix  D,  the  recursion  starts  with  L  =  1.6x64  +  40  =  142,  i.e.  ,  with 

—29 

J^^^(64)  =  0  and  (64)  =^10  .  The  maximum  term  in  the  power  series 

expansion  is  T^^  =  10  according  to  Table  4  in  Appendix  B,  so  no  loss 

of  accuracy  is  expected  from  this  effect.  However,  the  summation  of  the 
power  series  algorithm  ends  when  the  terms  become  smaller  than 
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10-^®  (=^30  approximately);  we  find  indeed  that  the  power-series  result 
is  accurate  to  30  decimal  places.  By  continuing  the  power-series  expan¬ 
sion  further  (to  60  terms),  the  more  accurate  result  of  the  footnote  on 
page  13  was  obtained. 
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Figure  1  -  Order-Argument  Domain  Covered  by  Algorithms  for  J^CX) 
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Figure  2  —  Results  Checked  or  Compared 


APPENDIX  A 


FORTRAN  IV  PROGRAM 

The  FORTRAN  IV  program  is  used  to  compute  Bessel  functions  of  the  first 
kind  of  positive  integer  order  and  positive  argument. 

The  program  is  written  as  a  double  precision  function,  with  indications 
of  the  changes  involved  when  a  single  precision  program  is  desired.  The 
program  has  been  checked  on  the  NSRDC  CDC  6700  system. 

The  calling  sequence,  variable  parameters,  and  options  are  explained 
in  comments  at  the  beginning  of  the  program. 
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C  PtSSEL  FUNCTION  ]Q3I,KOniJ 

C - - 

C 

DOUPLE  PRECISION  FUNCTION  DHSL  (LP » XP » I  OPT ) 

S  C 

C  BESSEL  FUNCTION  J(L,X)  OF  NON-NEGATIVE  INTEGER  OROER  L 
C  AND  NON-NEGATIVE  ARGUMENT  X. 

C 

C - 

10  c 

INTEGER  LP» IOPT»L» ISW»  TX, IM, I .K»KK»KM,LL  »LP1 ♦M 
REAL  ALM1,TX 

DOUBLE  RRECISION  XP . X » MACHFP *  PI » S ♦ S? » Y » Z » 77 
DOUBLE  PRECISION  A J . A JLM I ♦ A JL » A JLP 1 , AL » AL F , ALK , AM , AKK 
IS  C  INTEGER  lABS.  REAL  SORT,  DOUBLE  PRECISION  DABS»DCOS,DSIM.DSORT. 

C 

C - 

c 

C  ABSOLUTE  VALUES  OF  L  AND  X  ARE  TAKEN  BY  THIS  PROGRAM, 

20  X=DAPS(XP) 

L=IARS (LP) 

C 

C - 

C 

?S  C  MACHEP  IS  A  MACHINE  PRECISION  PARAMETER,  CHANGES  AFFECT  l.INES 
C  27»GHf72»  AND  86  IN  THIS  PROGRAM, 

MACHEP=2.<nM-97) 

C 

C - 

30  C 

C  THE  FOLLOWING  OPTIONS  ARE  AVAILABLE 
C  lOPT  =  1 

C  ALGORITHM  IS  BASED  ON  RECURRENCE  RELATION 

C  J(L-1»X)=(?L/X)*J(L»X)  -  J(L+1»X).  STARTING  WITH  SUFFICIENTLY 
3S  C  HIGH  L  FOR  WHICH  APPROXIMATELY  J ( L » X ) =MACHEP ♦  AND  J(L+1»X)=0, 

C  AND  NORMALIZED  BY  THE  SCALING  FACTOR  S=J  (0  »X) +2<»SUML  ( J  (2L  ♦X)  )  » 
C  WHICH  SHOULD  HAVE  EQUALED  ONE  FOR  THE  RIGHT  VALUES  OF  J. 

C  lOPT  =  2 

C  ALGORITHM  IS  BASED  ON  POWER  SERIES  EXPANSION 
40  C  SUMK((-1)»"K  (X/2)  (2K+L)  /  <  KFACT"  ( K +  L  )  F  ACT )  ) 

C  lOPT  =  3 

C  ALGORITHM  IS  BASED  ON  ASYMPTOTIC  SERIES  EXPANSION 
C  SORT(2/(PI<>X)  )  «(P(L»X)COS(ALF)  -  Q  (L  » X )  S I M  ( ALF )  ),  WHERE 
C  ALF  =  X  -  (.5L+.2S)PT»  AND,  ASYMPTOTICALLY,  WITH  MU=(2L)<^*2 
45  C  P(L»X)=1  -  (MU-1)  (MU-9)/{l*2  ^^(8X)^M^?)  +  .... 

C  Q(L»X)  =  {MU-1)/(1  ^^(BX))  -  (MU-1  )  (MU-9)  (MU-25)/(l<»2<»3  *(8X)^^^^3)) 
C  (-1)<>*INT(K/2)<^(MU-1)  ..  (MU-(2K-l)<»^^2)/(l<^2<^.,<^K  *(8X)<n^K)  +., 

C  lOPT  =  NOT  EQUAL  TO  1»2  OR  3 

C  ALGORITHM  DETERMINES  MOST  SUITABLF  OF  ABOVE  ALGORITHMS,  RUT 
50  C  MAY  LEAD  TO  UNSATISFACTORY  RESULTS,  ESPECIALLY  IF  L  AND  X 
C  ARE  OF  SAME  ORDER  AND  LARGER  THAN  100. 

C 
C 
C 

55  C 
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60 


C 

65  C 
C 
C 


C 

70  C 


C 

c 

c 

c 


80 


C 

C 

C 

C 


c 

c 

90  C 
C 


95 


100 


105 

C 


110 


PT=.314159?6535B9798?3B46^6433B33n+01 

AL  =  L 

AX  =  X 

IX  =  AX 

T5W=10PT 

TF  (TSW.LE.O)  ISW=4 
IF  (TSW.GF.4)  ISW=4 
GOTO  ( 10»?0»30,40) »  ISW 


40  FIMOS  MOST  5U1TABI..F  ALGORITHM  IF  POSSlHLF. 

40  IF  (IX.GF.50  .ANO.  4^^  I  X .  GF  .L<^L )  GOTO  30 
IMFOUAI.TTY  MAY  BF  PFOUCtn  TO 

IX,(.F.30  .ANO.  4■:^X.GF.L^n.  IK  M  ACHFP  =  ?^nM -47 )  . 

IF  (IX.I  T.l  .0-?.  ?<^TX.I..r.L)  GOTO  20 
IF  (TX.lT.lOO  .AMO.  L.LT.200)  GOTO  10 
TMFOUALITY  i^AY  BF  IMCPFASKD  TO 

IX.LT.I^B  .AMO.  I  .LT.200  IF  MACHFP  =  2*'if  (-47 )  . 

WRTTF {6» 15) 

SIGMAI.S  THAI  NO  SOITABLF  MFTHOO  IS  AVAILABLE^  RUT  PROVIDFS 
STIRLING  APPROXIMATION  TO  FIRST  TERM  OF  POWER  SERIFS. 

15  FORMAT (13H  WRONG  BESSEL) 

S=  (AX^^l  .  359  14091 42/ AL  )  ^^i^AL/SQRT  (6. 2831  8530  7<>-Al.) 

GOTO  47 


10  RECURPPMCF  RELATION. 

10  LPl=l.+  l 

KK=(B<UX)/5  +  40 

FINOS  ORDER  L  FOR  WHICH  APPROXIMATELY  J (L » X ) =MACHFP .  BFCOMES 
KK=(7^^IX)/5  +  25  IF  MACHFP=2^^<M -47 )  . 

K  =  2<MKK/?)  -1 

STARTS  NORMALIZING  SUM  S  (OF  TERMS  OF  EVFM  ORDER)  WITH 
HIGHEST  MON-ZERO  TERM. 

AJLP1=0. 

AJL=MACHEP 
S  =  Q, 

DO  23  1=1 »K»2 

LL=K-I+2 

AL=LL 

AJLMl  =AJI.■!^2.^l■AL/X  -  AJLPl 

AJLP1=AJL 

AJL=AJLM1 

TF  (LL.FO.LPl)  AJ=AJL 

LL=K-I+1 

AL=LL 

AJLM1  =  AJL<^2.<^AL/X  -  AJLPl 
S=S+AJLP1 

ADDS  THE  NEXT  LOWER  FVFN  ORDER  lERM  TO  THE  NORMALIZING  SUM  S. 
AJLP1=AJL 
AJL  =  AJI..M1 

IF  (LL-.EO.LPl)  AJ  =  AJL 
23  CONTTMUF 
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c 


1  15  c 
C 
C 

C  20 
20 

120 


33 

125  C 

37 


130  57 


135 

C 

140  C 
C 

C  30 
30 

145 

C 

150 

77 

155 

160 


S-?.<^S  +  AJL 

COMPLETFS  the:  WORMALIZIMF,  sum  S=J(0»X) +2^^SUML(J(2L^X)  )  . 

S=AJ/S 

GOTO  47 


PnWFP  SFPTETS  EXPANSION. 

Y=X/2. 

ALF=1  . 

IF  (L.FO.O)  GOTO  37 
00  33  1M=1»L 
AM=IM 

ALF=ALF«AM 

ALF  IS  (K+L) -factorial.  NOW  SET  UP  FOR  AND  EXECUTE  SUMMATION. 
Z=1 ./ALF 

S=7 
KK  =  1 
AKK=KK 
ALK=L+KK 

7  =  -z^>Y^^Y/(AKK^^ALK) 

S  =  S  +  7 
ZZ  =  f)ABS(Z) 

IF  (ZZ.LT.MACHEP)  GOTO  47 
IF  (KK.GF.60)  GOTO  47 
KK  =  KK+  I 
GOTO  57 


ASYMPTOTIC  SERIES  EXPANSION. 

Y  =  B.<^X 

ALM1=4<:l<^L-1 
TX  =  ?.«X  +  l  . 

KM=.5<MTX  +  SQRT  (TX-^TX  +  ALMl  )  ) 

FINOS  MINIMUM  TERM  IN  CASE  IT  OCCURS  BEFORF  MACHEP  IS  REACE^F^, 
AM=4<H.-‘^L 
Z  =  l. 

S=l. 

S2  =  0. 

KK=1 

AKK=KK 

M  =  4<^KK<:-KK  -  4<!-KK  +  1 
ALK=M 

Z=Z<MAM-ALK)  /  (AKK^^Y) 

S2=S2+Z 

IF  (KK+l.GF.KM)  GOTO  167 
akk=kk+i 

ALK=M  +  e^^KK 
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APPENDIX  B 
POWER  SERIES  TABLES 
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TABLE  2  -  POWER  SERIES  FINAL  TEEM  AND  COMMON  LOGARITHM  OF  ITS  MAGNITUDE  FOR  MACHINE  PRECISION  10 

(Algorithm  continues  for  60  terms  or  until  required  precision  is  reached) 
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TABLE  4-POWER  SERIES  MAXIMUM  TERM  AND  COMMON  LOGARITHM  OF  ITS  MAGNITUDE 

(Dashed  line  indicates  boundary  at  which  loss  of  accuracy  may  be  2  decimals; 
algorithm  chooses  power  series  if  X  is  less  than  ^  L  or  if  X  is  less  than  1^ 
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TABLE  5-ASYMPTOTIC  SERIES  MAXIMUM  TERM  AND  COMMON  LOGARITHM  OF  ITS  MAGNITUDE 

(Dashed  line  indicates  boundary  at  which  loss  of  accuracy  may  be  2  decimals; 
algorithm  chooses  asymptotic  series  if  X  is  greater  than  L^/4) 


<N 

O 


CNl 

iH 


vO 

m 

CN| 


00 

rH 


<f 

VD 


CO 


vO 

rH 


00 


ix^ 


m 


CO 

I 


iH  CN 


CO 

I 


rH 


CM 

I 


rH  VO 


CM 

I 


rH  CO 


CM 

I 


rH  O 


CM 

I 


rH 

I 


H  -d- 

rH 

I 


rH  rH 

iH 

I 


00 

I 


H  to 

r 


rH 


CM 

I 


rH 


CM 

I 


CM 

I 


00 

• 

?— I 
I 


rH 

( 


« 

I 


CM 

I 


00 
I— ! 

I 


to 

rH 

I 


CM 

rH 

I 


rH  0^ 


vO 

I* 


rH  to 


rH  CM 

rH 

I 


rH  a\ 


VO 

I 


CO 

I* 


rH  O 


Ov 

r 


vO 

I 


CO 

I 


H  CO 
I 

H  O 


H  CO 


CO  O 


rH  O 


rH  CO 


CO  O 
rH 


VO 

CM 


J 


■H  fO 


r 


J 


in  CO 
rH  • 

in 


CO  O 


CO 

r 


rH  CO 


r; 


J 


vO 

CM 


rv  vO 
CM 


CO  m 

rH  • 

in 


O  sd- 
CM  • 

o 


VO  CJ\ 
CM  • 


ov 

CO  • 
rH 
CM 


O  CO 

m  • 

in 

CO 


VO 

CM 


to  ov 
rH  • 

m 


o  tn 

CO  • 
CM 


CO 

in  • 

CM 


C3V  tn 


o  vo 

o  • 


CO  00 


CO 

o 


rH 

CO  • 
CM 


o  o 

vO  • 
vO 
CM 


VO  m 

o  • 

iH  o 
m 


00  CO 

m  • 

H  o 
a\ 


o  m 

o  • 

CM  <t 


vo  O 
CM  • 
CM  a\ 

o 

CM 


O 

^  • 

CM  a\ 

CM 


O  H 
CM 


O 
CM 


tH 

CO 

^  vO 

m 

o 

CM 

CM 

CO 

CO  <1* 

rH 

• 

CO  • 

CM 

CM 

vo  tn 

O 

vo 

rH 

CO 

vO 

OV 

a\  in 

rH 

• 

<X\  • 

CO 

rH 

CM 

00 

00 

rH 

m 

O 

tn 

CO 

O 

• 

o  • 

o 

o\  o 

<JV 

CM 

00 

CM 

CM  O 

in 

• 

VO  • 

Ov 

0\  CM 

rH 

CM 

rH 

rH 

rH 

CM 

CM 

CO 

• 

o\  • 

o 

CJn  vo 

vO 

rH 

m 

'd- 

rH 

vo 

00  H 

OV 

• 

O  • 

O  00 

O 


H  H 
rs. 


1  \o 

iH  O 

CO 

O 

O  CO 

m  CM 

VO 

O  O 

00  rH 

CM 

VO  OV 

• 

• 

rH  • 

CM  • 

in  • 

CM  • 

SJ*  • 

o  • 

H  • 

1 

rH  ’ 

rH 

rH  C7V 

CM  CO 

in  00 

O  (N 

1 

rH 

m 

CO 

tn 

in 

H  CM 

rH 

CO 

00 

O 

J 

CM 

CO 

iH  CO 

1 

i-H 

CM  CM 

00  CM 

o  o 

-vt"  o^ 

CM  <1“ 

00  vo 

o  <■ 

• 

• 

♦ 

rH  • 

CM  • 

vo  • 

CM  • 

m  • 

o  • 

CM  • 

1 

CM 

00 

in 

fH  in 

CM  00 

in  o 

o  o\ 

CM 

vO 

CM 

rH 

H  CM 

rH 

o 

CO 

rH 

CM 

o 

CM  00 

vD 

00 

CM 

o  o 

CM  <1- 

vo  in 

vO 

O  Oi 

CM  00 

• 

• 

♦ 

rH  • 

CO  • 

VO  • 

CM 

tn  * 

rH  • 

CM  • 

1 

CO 

CM 

rH  CO 

CM  <1* 

m  CO 

O  vO 

rH 

CO 

00 

iH 

o 

vO 

rH  CO 

CM 

in 

rH 

vO 

rH 

CM 

CO 

CO  vo 

00 

in  vo 

rH  CM 

CO  CM 

r-  vo 

in  CM 

rH  vO 

CO  vO 

• 

• 

• 

rH  • 

CO  • 

VO  • 

CM  • 

in  • 

rH  • 

CM  • 

rH 

in 

vD 

CO 

VO 

rH  rH 

CM  rH 

in 

o  <• 

rH 

sT 

O 

in 

00 

rH 

rH  -d- 

iH 

CM 

in 

CO 

CJV 

rH 

CM 

CM 

00 

vO 

CM 

00 

vo 

CM 

<1- 

rH 

CO 

vD 

CM 

m 

rH 

CM 

rH 

CM 

tn 

o 

rH 

28 


TABLE  6-ASYMPTOTIC  SERIES  MINIMUM  TERM  AND  COMMON  LOGARITHM  OF  ITS  MAGNITUDE 

(Dashed  lines  indicate  boundaries  at  which  magnitude  attains  single  machine  precision  (10 
and  double  machine  precision  (10  algorithm  chooses  asymptotic  series  if  X  is  greater 
than  30— single  precision— or  if  X  is  greater  than  50— double  precision) 
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TABLE  8-ASYMPTOTIC  SERIES  FINAL  TERM  AND  COMMON  LOGARITHM  OF  ITS  MAGNITUDE  FOR  MACHINE  PRECISION  10 
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APPENDIX  D 


RECURRENCE  RELATION  TABLE 

TABLE  9-BESSEL  FUNCTIONS  WHOSE  MAGNITUDE  APPROXIMATELY  EQUALS 

MACHINE  PRECISION 

(Fast  rule:  L  w  1.4X  +  25,  mach.  prec.  10  L  «  1.6X  +  40, 
mach.  prec.  10“^®) 
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Log  Jl(X) 
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Log  Jl(X) 
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10 
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-28.3 

20 
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20 
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30 
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40 

103 

-28.6 
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50 

119 

-29.0 

60 

109 

-13.7 

60 

134 

-28.9 

70 

123 

-13.7 

70 

149 

-29.0 

80 

137 

-13.8 

80 

163 

-28.7 

90 

151 

-13.8 

90 

178 

-29.0 

100 

165 

-13.9 

100 

192 

-28.8 

200 

302 

-13.8 

200 

332 

-28.8 

400 
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